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Abstract The rapid racerunner, Eremias velox, is a widely distributed lizard from the northern Caucasus across entire 
Central Asia eastward to China. It is increasingly common to accept E. velox as a species complex in its entire range. 
To date, published morphological and molecular systematic hypotheses of this complex are only partially congruent, 
and its taxonomic status and evolutionary history are still far from clear. The mitochondrial cytochrome b gene and 
12S rRNA sequences were used to evaluate the taxonomy of this complex, with particular attention to the phylogenetic 
placement of populations in northwestern China. Examination of the phylogenetic analyses recovers seven distinct, 
biogeographically discrete, and well-supported clades, revealing genetically identifiable populations corresponding to 
some previously morphology-defined subspecies. Chinese E. v. roborowskii appears to have split from other Central 
Asian rapid racerunner lizards well before differentiation occurred among the latter taxa. Specifically, we corroborate 
that there are two subspecies occurring in China, i.e., E. v. velox and E. v. roborowskii. We recommend a novel 
subspecific status for the phenotypically and genetically distinct populations in southern Aral Sea region of Uzbekistan 
previously assigned to Æ. v. velox. Finally, each of the three independently evolving lineages from Iranian Plateau should 
be recognized as three species new to science under the general lineage concept. 


Keywords mtDNA, subspecies, Eremias velox complex, taxonomy, phylogeography 


1. Introduction been assigned to the genus, and much of their phylogeny 
and taxonomy is controversial (e.g., Szczerbak, 1974; 
Racerunner lizards of the genus Eremias (family 
Lacertidae) are the dominant reptiles in the deserts and 
steppes of Central Asia. The taxonomy of this genus has 
a long complex history, and the delimitation of species/ 
subspecies of Eremias is known for being notoriously 
difficult. This is, in part, caused by the great geographic 
variation in most morphological characters in many of 
the species (e.g., Szczerbak, 1974, 1975; Ananjeva et al., 
1998; Chirikova, 2004). Currently, over 36 species have 


Guo et al., 2011). Among the many controversies, the 
taxonomy of the Eremias velox complex is perhaps one 
of the most confusing. The distribution of this group 
ranges from the northern Caucasus across entire Central 
Asia eastward to China (Figure 1). Traditionally, the rapid 
racerunner was considered as a single species with three 
or four described subspecies (Szczerbak, 1974, 1975; 
Eremchenko and Panfilov, 1999). However, as noted by 
Rastegar-Pouyani (2009), it is increasingly common to 
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accept E. velox as a clade in its entire range representing a 
species complex. The nominate form occurs in the larger 
parts of the distribution range. The Caspian subspecies Æ. 
v. caucasia occurs at the western part of the range. The 
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subspecies E. v. roborowskii is endemic in Turpan-Hami 
Depression and Dunhuang Basin in northwestern China. 
Recently Eremchenko and Panfilov (1999) described the 
fourth subspecies, E. v. borkini from the highlands of 
Issyk-Kul lake depression in Kyrgyzstan. 

Many recent taxonomic treatments reveal varied 
opinions as to the number of species. Some herpetologists 
recognized a single widespread species (e.g., Szczerbak, 
2003; Ananjeva et al., 2006; Sindaco and Jeremcenko, 
2008), whereas others (Rastegar-Pouyani, 2009; Rastegar- 
Pouyani et al., 2012) hypothesized the existence of at least 
four distinct species and three well-supported subspecies 
(Figure 2). Some of the differences of opinion result 
from the accumulation of mitochondrial DNA sequence 
information (Rastegar-Pouyani ef al., 2012). At present 
it seems reasonable that mtDNA sequence information 
supports the recognition of the Iranian populations of 
E. velox as three distinct species based on differences in 
distribution, morphology, habitat preference, and mtDNA 
sequence divergence (Figure 2). 

The subspecific status of E. v. roborowskii is still 
controversial, although Szczerbak (1975) confirmed 
its status according to color, pattern and characters 
of pholidosis. Boulenger (1921) suggested that the 
morphological variation in the populations of E. v. 
roborowskii 1s not sufficient to support its subspecies 
status. The validity of this subspecies has not been 
recognized for a long time (Mertens and Wermuth, 1960). 
More recently, on the basis of multivariate morphological 
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analyses, Wang et al. (2014) argued that there is no 
subspecies differentiation between populations from 
Turpan-Hami Depression and from Dzungarian Basin 
and Ily Valley in China. Thus, Wang et al. (2014) implied 
that E. v. roborowskii may be a synonym of nominate 
subspecies. 

Establishing correct species/subspecies limits depends 
on adequate sampling of populations throughout the 
range. Given the previous studies, it might appear that 
taxonomic limits were clear. Rastegar-Pouyani (2009) 
made the first attempt to elucidate the phylogenetic 
relationships among the rapid racerunner in Iran and 
Central Asia on the basis of ISSR-PCR fingerprints 
data. The results were tentative due in part to limited 
taxonomic sampling. As shown by simulations and 
empirical studies, more taxa and data do affect the 
phylogenetic reconstruction (e.g., Pollock et al., 2002; 
Zwickl and Hillis, 2002). Rastegar-Pouyani et al. (2012) 
studied patterns of mitochondrial DNA cytochrome b 
(cyt b) and 12S rRNA variation in 70 specimens from 
13 geographically distant localities in Iran and Central 
Asia. Their results contradicted previous systematic 
hypotheses; revealed surprising relationships between 
Iranian and Central Asiaian lineages and uncovered novel, 
cryptic evolutionary lineages (i.e. new putative species). 
Moreover, they argued that the eastern and southern 
Kazakhstan populations correspond to the traditional Æ. 
v. roborowskii. However, Rastegar-Pouyani et al. (2012), 
who used most samples already analyzed by Rastegar- 


Figure 1 Entire distribution range of Eremias velox complex showing traditional subspecies’ range with different background color, and the 
sampling localities for different populations (numerical labels 14-30 with green dots), white dots with numbers 1-13 representing those from 
Rastegar-Pouyani et al.(2012). The locality number corresponds with that in Appendix A. 
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Figure 2 Rastegar-Pouyani et al.’s (2012) phylogenetic hypothesis for populations of the Eremias velox complex from Iranian Plateau 
and Central Asia using partitioned Bayesian analyses of mtDNA data. A, B, C, D and E indicate the major clades. Values beside the node 


represent Bayesian posterior probabilities. 


Pouyani (2009), were unable to sample adequately the 
eastern part of the taxon’s range. There was no sample 
from China included in their work. In addition, Rastegar- 
Pouyani et al. (2012) speculated that E. v roborowskii 
occurs in western China, eastern Kazakhstan and Eastern 
Uzbekistan. This inappropriate interpretation may give 
confusing clues about the identity of this subspecies, 
resulting in incorrect inference of its phylogenetic 
placement. It appears likely that the inclusion of more 
rapid racerunner samples may get a better illustration of 
their phylogeographic structure, and may also provide 
valuable information for their taxonomy. 

Thus, as an extension of the molecular phylogenetics 
of the rapid racerunner, an effort has been made to collect 
some specimens from northwestern China and Russia. 
We further determined the complete cyt b gene and 12S 


rRNA segment to gain a more comprehensive view on 
phylogeographic patterns in Æ. velox complex. This study 
specifically aims to: (1) re-evaluate the phylogeographic 
structure of E. velox complex, (2) test the phylogenetic 
affinities of the rapid racerunner populations in China, 
and (3) test the validity of the subspecific status of E. v. 
roborowskii. The hypothesis that E. v. roborowskii is valid 
particularly predicts that it will correspond with distinct, 
well-supported, and moderately divergent mtDNA 
haplotype clade. 


2. Materials and Methods 


2.1 Taxon sampling We collected a total of 87 samples 
representing natural populations at 17 localities. These 
included 57 samples from 12 localities in China, 30 
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samples from five localities in Russia, covering a 
large territory occupied by E. velox in the two regions 
(Appendix A, Figure 1). To ensure accuracy more than 
four representatives of each population were collected. 
However, only a single individual was collected from 
six sites due in part to sampling difficulties. Populations 
were assigned subspecific status based on location and 
phenotype. Due to logistical and legal obstacles to filed 
work, no sample of E. velox borkini from the highlands of 
Kyrgyzstan was available. Based on present knowledge 
of phylogenetic relationships among Eremias lizards, E. 
argus and E. persica were chosen as outgroup taxa (Guo 
et al., 2011). All voucher specimens were deposited in 
the herpetological collections of Chengdu Institute of 
Biology, Chinese Academy of Sciences or the Department 
of Zoology, University of Guelph. Detailed information 
about all specimens and sequences used in this study is 
listed in Appendix A. 


2.2 DNA extraction, amplification, and sequencing 
Genomic DNA was extracted from liver or tail muscle 
tissues stored in 95% ethanol following the protocol of 
Aljanabi and Martinez (1997). Complete sequence of 
the mitochondrial DNA cyt b gene was amplified using 
the primers L14919 and H16064 from Burbrink et al. 
(2000). A ~371bp fragment of 12S rRNA was amplified 
using the primers L1091 and H1478 from Kocher et al. 
(1989). Each of our 50 ul PCR reactions contained 25 ul 
of 2 x EasyTaq SuperMix, 0.4 uM of each primer, and 
1-2 ul genomic DNA. The PCR protocol involved initial 
denaturation at 94°C for 180 s followed by 30 cycles of 
94°C for 30 s, 43°C for 30 s (50-52°C for 12S rRNA), 
and elongation at 72°C for 90 s (60 s for 12S rRNA); 
and final extension at 72°C for 8 min. Negative controls 
were run for all amplifications. Amplified products were 
visualized on 0.8% agarose gels and purified using the 
DNA gel extraction kit (Omega Bio-Tek). The purified 
products were sequenced directly by Majorbio or Sangon 
(Shanghai) Co., Ltd., using the same primers as PCR for 
12S rRNA, and L14919 from Burbrink et al. (2000) in 
combination with L14761 as well as H14892 from Zhao 
et al. (2011) for cyt b. All novel sequences are deposited 
in GenBank with accession numbers KF999318- 
KF999491. 


2.3 Sequence alignment and analyses We have 
retrieved sequences of cyt b and 12S rRNA of 39 
individuals from GenBank, which included sequences of 
37 individuals published in Rastegar-Pouyani et al. (2012) 
and sequences of two individuals used as outgroup in this 
study (Appendix A). Sequences of cyt b were translated 


to amino acid sequences using MEGAS (Tamura et al., 
2011) to verify the data. A total of 126 sequences of 
cyt b and 12S rRNA were first aligned using Clustal X 
2.0 (Larkin et al., 2007), respectively, with default gap 
penalties. Then the aligned matrice were re-checked 
visually. 

Phylogenetic congruence of cyt b and 12S rRNA data 
sets were tested by the partition homogeneity test of Farris 
et al. (1995) with PAUP* 4.0b10 (Swofford, 2003). The 
partition homogeneity test supported the combination of 
the cyt b and 12S rRNA data sets (P = 0.1). Homogeneity 
of base frequencies was evaluated using Chi-square °) 
tests implemented in PAUP*. Base composition at the 
third position exhibited no heterogeneity: the combined 
data, y = 16.66, df = 237, P = 1.00; first position, y = 
7.18, df = 237, P = 1.00; second position, y’ = 2.02, df 
= 237, P = 1.00; third position, y’ = 48.99, df = 237, P 
= 1.00; and 12S rRNA, 7’ = 12.11, df = 237, P = 1.00. 
Given these results, a sensible approach to analyze these 
data phylogenetically would be to apply a time-reversible 
Markov model. 


2.4 Phylogenetic analyses To estimate phylogenetic 
relationships, Bayesian inference (BI), maximum 
likelihood (ML), and maximum parsimony (MP) were 
implemented for the combined data. Each haplotype 
was treated as a taxon in the analyses. MP analyses with 
PAUPRat (Sikes and Lewis, 2001), based on Parsimony 
Ratchet (Nixon, 1999), were conducted using 1000 
ratchets with 10 iterations per replicate. The shortest 
equally most parsimonious trees were combined to 
produce a Strict consensus tree. Consistency (CI), 
retention (RJ), and rescaled consistency (RC) indices were 
used to evaluate homoplasy in character data and overall 
support for optimal trees. Support for individual branches 
was evaluated using non-parametric bootstrapping 
(Felsenstein, 1985) in PAUP* with 1000 pseudoreplicates, 
TBR branch swapping, simple taxon addition with one 
tree held at each step, and a maximum of 100 trees saved 
per replicate in order to decrease the time needed to run 
large bootstrap replicates. 

To infer relationships under BI we used MrBayes 
v3.2.2 (Ronquist et al., 2012). We optimized several 
parameters prior to conducting our final analyses using 
the Markov Chain Monte Carlo (MCMC) method 
implemented in MrBayes, including partitioning strategy, 
models of molecular evolution, length of MCMC 
analyses. We assessed two alternative partitioning 
strategies: (1) a single partition and (2) a four-partition 
strategy that applied a single partition to each of the three 
codon positions of cyt b and a fourth partition for 12S 
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rRNA sequence. 

The Akaike Information Criterion (AIC; Akaike, 
1974), as implemented in MrModeltest v2.3 (Nylander, 
2008), was used to select the appropriate model of 
sequence evolution for each partition. To compare the 
two alternative partitioning strategies, we conducted 
two independent analyses of each strategy in MrBayes. 
Each analysis consisted of three heated chains and one 
cold chain, with an incremental heating temperature 
of 0.02, and an exponential distribution with a rate 
parameter of 25 as the prior on branch lengths (Marshall, 
2010). All analyses were run for 20 million generations 
with parameters and topologies sampled every 1000 
generations. To ensure MCMC convergence, a burn-in 
period for all analyses was diagnosed in two ways: (1) 
by the average standard deviation of split frequencies 
(ASDSF; Lakner et al., 2008) between two MCMC 
analyses run independently, with levels below 0.01 
being considered indicative of convergence, and (2) by 
the potential scale reduction factor (PSRF; Gelman and 
Rubin, 1992) comparing the variance within and between 
runs, with levels approaching 1.0 as runs converge. To 
compare performance of the two partitioning strategies 
we calculated Bayes factors from the marginal likelihood 
scores estimated by the stepping-stone method (Xie et al., 
2011) in MrBayes 3.2.2. Comparison of alternative 
partitioning strategies strongly favored four partitions 
over one partitioning scheme (2Ln Bayes factor > 600). 
After selecting the optimal partitioning strategy, we 
conducted two analyses with different starting trees 
using the optimized parameters. We then used the sumt 
command in MrBayes to generate a consensus topology 
and associated Bayesian posterior probabilities (BPP) 
from a pooled post-burnin sampled from two analyses. 

Partitioned maximum likelihood (ML) analyses were 
conducted in RAxMLHPC v7.0 (Stamatakis, 2006) on 
the combined dataset with the same partitioning strategy 
as for Bayesian inference. The more complex model 
(GTR +/+) was used for all subsets, and 100 replicate 
ML inferences were performed for each analysis. Each 
inference was initiated with a random starting tree 
and nodal support was assessed with 1000 bootstrap 
pseudoreplicates (Stamatakis et al., 2008). 


2.5 Bayesian tests of topology hypotheses Bayes 
factors were used to compare the unconstrained 
Bayesian tree topology to Bayesian trees with ‘hard’ 
constraints. Here we tested the following taxonomy 
and geography-based hypotheses to evaluate the 
taxonomic status within E. velox complex: (1) 
Are populations within the historically accepted 


distributions in northwestern China monophyletic? 
(2) Are the Iranian populations monophyletic? (3) Do 
Iranian populations belong to nominate subspecies? 
(4) Are populations of E. v. roborowskii and nominate 
subspecies monophyletic? We constrained the topology 
for each of the above putative groupings and performed 
a stepping stone run of 10 million generations (50 
steps with convergence being reached in each step) 
from which we obtained a mean marginal likelihood 
estimate. Convergence was examined through diagnostic 
plots of standard deviation of the split frequencies 
and screening similarity in the two independent 
marginal likelihood estimates (Ronquist et al., 
2011). Each of these estimates were compared using 
Bayes factors (Kass and Raftery, 1995) to the marginal 
likelihood estimate obtained from the unconstrained 
MrBayes run of the same length. The stepping stone 
function (Xie et al., 2011) implemented in MrBayes 
v3.2.2, offers an improved estimation of marginal 
likelihood over harmonic mean estimation. 


3. Results 


3.1 Sequence characteristics Two mitochondrial genes 
(cyt b and 12S rRNA) were obtained, yielding alignments 
of 1143 and 372 bp, respectively. 126 sequences 
(including two outgroup taxa) revealed 80 haplotypes. 
Alignment was unambiguous and 1 indel was inferred 
in 12S rRNA between outgroup and ingroup taxa. No 
premature stop codons were observed in cyt b, indicating 
that the obtained sequences were mitochondrial in origin 
and not nuclear pseudocopies. Average base composition 
is A: 27.55%, C: 30.69, G: 14.72%, T: 27.04%, when the 
outgroup taxa are combined. Nucleotide composition 
showed an anti-G bias and A+7-richness (54.59%), which 
is a characteristic of the mitochondrial genome. Variable 
and parsimony-informative characters are: 441 and 319 of 
1143 (cyt b); 84 and 57 of 371 (12S rRNA). 

Ax test at the 5% level of significance for differences 
in base frequencies showed that there was no base 
compositional heterogeneity among sequences, which is 
known to adversely affect phylogenetic inference (Jermiin 
et al., 2004). 


3.2 Phylogenetic analyses The heuristic search of the 
combined data resulted in 914 equally parsimonious trees 
of 1136 steps, with moderate CI (0.579) and high RI 
(0.887). These are presented as a majority-rule consensus 
tree where branches with bootstrap support lower than 
50% are collapsed (data not shown). AIC analyses 
conducted with the aid of MrModeltest identified the 
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HKY + Z + T model as the most appropriate model for 
partition of cyt b 1" codon position, and the GTR + J + 
T model for partitions of cyt b 2" codon position, cyt b 
3" codon position, 12S rRNA or combined data. Bayesian 
analyses were conducted for 20 million generations, with 
a conservative burn-in of six million generations that 
provided dense sampling of the posterior distribution. The 
MP, ML and Bayesian trees differed in support values, 
but not in the compositions of the major lineages. The 
bootstrap support was also marked on the branches that 
receive such support in the Bayesian tree (Figure 3). 
The phylogenetic trees recovered the monophyly of the 
E. velox complex with high support values. As shown 
in Figure 3, seven major clades connected deep in the 
complex phylogenetic history were identified in all 
reconstructions with strong posterior probabilities and 
bootstrap values, which generally correspond well to 
geographical regions within the distribution range. These 
clades were designated I, II, III, IV, V, VI and VII in 
Figure 3. The subdivision and phylogenetic positions of 
the Iran populations (J, II, III) are completely identical 
with the clades A, B, C of Rastegar-Pouyani et al. 
(2012), respectively. The basal split in the complex 
is between the northern Iranian population (clade I) 
and a highly supported clade of the other populations 
from Iran, Central Asia, China and Russia. However, 
clade IV, corresponding to E. v. roborowskii, covers the 
populations from Turpan-Hami Depression, and forms 
the sister taxon to the other populations excluding those 
from Iran. The clade V, corresponding to EF. v. caucasia, 
covers populations from North-eastern Russian Caucasus 
and Western Kazakhstan. The clade VI is composed of 
populations from Uzbekistan, conventionally assigned to 
nominate subspecies. The clade VII covers populations 
from Eastern Kazakhstan and Dzungarian Basin as 
well as Ily Valley in China, which corresponds to 
nominate subspecies. Overall, the notable differences 
between Rastegar-Pouyani et al. (2012) and our results 
presented here lie in the phylogenetic placement of E. v. 


roborowskii. 


3.3 Bayesian tests of topology hypotheses We compared 
the marginal likelihood estimates of topologically 
constrained and unconstrained run of MrBayes for 
possible monophyletic group within E. velox complex 
(Table 1). In all cases, there was very strong (Bayes 
factor > 100) evidence against the constrained topologies. 
Thus, four alternative phylogenetic hypotheses were 
significantly rejected, such as a monophyletic group of 
Chinese rapid racerunner populations, monophyly of 
Iranian populations, assignment of Iranian populations to 
nominate subspecies, and assignment of E. v. roborowskii 
to nominate subspecies. 


4. Discussion 


This study provided strong support for the monophyly of 
E. velox, consistent with the results of Guo et al. (2011) 
and Rastegar-Pouyani et al. (2012), and thus providing a 
firm basis for future testing of alternative biogeographic 
models. The Iranian populations were recovered here as 
a paraphyletic assemblage, representing three distinct, 
biogeographically discrete, well-supported clades that are 
not each other’s closet relatives. This is in accord with the 
results of Rastegar-Pouyani et al. (2012). 

The morphological differences (Rastegar-Pouyani 
et al., 2012), allopatric ranges and substantial genetic 
divergence in the mtDNA bear evidence of long standing 
isolation between several populations. The haplotype 
clades associated with each lineage on Iranian Plateau 
are differentiated by mean uncorrected levels of sequence 
divergence that often exceed 10% (Table 2). This degree 
of divergence suggests that evolutionary separation of 
the three clades occurred millions of years ago, likely 
sometime during the Miocene based on a molecular clock 
calibration of 1.25% pairwise divergence per million 
years for mtDNA within the Æ. velox complex (Rastegar- 
Pouyani et al., 2012). Concordance between strongly 
differentiated mtDNA haplotype clades and phenotypic 


Table 1 Bayes factors of stepping-stone-based estimates of marginal likelihood for comparisons of alternative phylogenetic hypotheses. 
2In Bayes factors >10 are considered very strongly different (Kass and Raftery, 1995), indicating evidence against alternative hypotheses. 


. . a Ln Marginal likelihood 

Alternative phylogenetic hypotheses Description of constraint 2Ln Bayes factor 

LnZ: unconstrained LnL: constrained 

H1 pre Ea oE fápid radenintet populationsin. a China -8022.39 -8121.34 197.9 

Ww Monophyly of rapid racerunner populations in All samples from Iran: -8022.39 8076.55 108.32 
Iran clade I, II, III 

H3 Nominate subspecies affinity of Iranian rapid clade I, Il, IL, VI, VII -8022.39 -8082.49 120.2 
racerunner populations 

H4 Monophyly of E. v. velox and E. v. roborowskii clade IV, VI, VII —8022.39 -8090.19 135.6 
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Figure 3 Hypothesized relationships of Eremias velox complex included in this study, illustrated by the majority-rule consensus tree 
resulting from partitioned Bayesian analyses. Bayesian posterior probabilities, maximum parsimony and likelihood bootstrap values are 
shown. Nodes supported by > 0.95 Bayesian PP and > 70% BP were considered highly supported. Terminals are labeled with haplotypes. 
Roman alphabet labels correspond to clades referred to the Results and Discussion. Representative photographs of female adults show the 
dorsal and lateral differentiation among three subspecies. 
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Table 2 Genetic distances within and between outgroup species and FE. velox complex included in this study. Mean uncorrected (p) distances 
are below diagonal; distances calculated using MEGA’s Maximum Composite Likelihood model with rate variation among sites defined by 
the gamma shape parameter (0.2366) are above the diagonal. Mean maximum likelihood divergences within each species or subspecies are in 
the diagonal with bold text. 


1 2 3 4 5 6 7 8 9 
1 E.argus NA 0.419 0.483 0.442 0.443 0.454 0.453 0.471 0.433 
2 E. persica 0.161 NA 0.416 0.442 0.489 0.462 0.441 0.433 0.401 
3 clade I 0.172 0.158 NA 0.173 0.251 0.267 0.272 0.325 0.299 
4 clade II 0.164 0.159 0.099 0.03 0.092 0.145 0.152 0.171 0.164 
5 clade III 0.165 0.17 0.122 0.066 0.033 0.14 0.167 0.186 0.17 
6 clade IV 0.17 0.165 0.128 0.092 0.089 0.003 0.152 0.17 0.153 
7 clade V 0.168 0.161 0.129 0.095 0.101 0.094 0.002 0.111 0.093 
8 clade VI 0.166 0.154 0.14 0.099 0.105 0.098 0.075 0.006 0.04 
9 clade VII 0.162 0.152 0.135 0.097 0.1 0.093 0.067 0.033 0.012 


variation (Rastegar-Pouyani et al., 2012) supports the 
hypothesis that each of the three independently evolving 
lineages from Iranian Plateau deserves recognition at the 
species-level under the general lineage concept of species 
(de Queiroz, 1999; 2007) and associated operational 
species-delimitation criteria (Wiens and Penkrot, 2002). 
Rastegar-Pouyani et al. (2012) demonstrated the 
basal positions of the Iranian clades for the first time 
and considered them three distinct species. However, 
their inference of the phylogenetic placement of E. v. 
roborowskii was misleading. In support of Rastegar- 
Pouyani et al. (2012), our results confirm the basal 
position of the Iranian clades and statistically reject the 
placement of Iranian clades into the nominate subspecies 
from Central Asia. Surprisingly, we find that E. v. 
roborowskii is the sister taxon to E. velox caucasia plus 
E. v. velox, i.e., other populations from Central Asia. 
Our results clearly confirm that all of the major clades 
inferred from mitochondrial haplotypes are concordant 
with geographical regions. Iranian clades at the basal 
position of the phylogenetic tree may imply that E. velox 
originated from the Iranian Plateau and invaded Central 
Asia as Rastegar-Pouyani et al. (2012) suggested. If we 
take this scenario as a reliable explanation for the origin 
and evolution of E. velox, we consider the clade IV (E. 
v. roborowskii) as first split from entire Central Asian 
populations, resulting from the barrier effect of Tianshan 
Mountains with their dramatic uplift since Tertiary (Macey 
et al., 1999; Guo and Wang, 2007; Zhang et al., 2008; 
Melville et al., 2009; Solovyeva, 2013). Geographically, 
with the uplift of Tianshan Mountains, Turpan-Hami 
Depression and Dzungarian Basin were blocked from 
each other, resulting in a historical vicariant pattern and 
allopatric distribution between E. v. roborowskii and E. v. 
velox as well as E. v. caucasia. Wang et al. (2014) failed 
to recognize the validity of E. v. roborowskii due to a 


lack of understanding the discrete characters between E. 
v. velox and E. v. roborowskii, including dorsal pattern 
and lateral color pattern. In general, the adults of Æ. v. 
roborowskii have a dorsal pattern of irregular dark spots, 
with rows of bright blue-edged, dark eye-like spots on the 
lateral sides (Zhao, 1999; Szczerbak, 2003). 

In addition, clade V, corresponding to traditional Æ. 
velox caucasia, forms the sister taxon to Central Asian 
clades, in accord with the results of Dolotovskaya et 
al.(2007) inferred from ~ 1800 bp mtDNA of cyt b, COI 
and 16S rRNA . Rastegar-Pouyani et al. (2012) speculated 
that a single specimen from Western Kazakhstan 
belongs to E. velox caucasia. The affinity of that single 
specimen was corroborated in this study for there to be 
strong support for its nesting inside clade V. The Central 
Asian clades, corresponding to conventional nominate 
subspecies, comprise clades VI and VII. The patterns of 
mtDNA variation support a sister relationship between 
clade VI and clade VII. Clade VI comprises samples from 
Uzbekistan, whereas clade VII consists of populations 
from Kazakhstan and Dzungarian Basin and Ily River 
Valley, Xinjiang, China. Judging from the type locality 
of the nominate form of E. velox, Inderskija mountain 
in Kazakhstan, we consider clade VII as nominate 
subspecies. Szczerbak (1975) examined the relationships 
of E. v. roborowskii from Ssatschsheu (Dunhuang Basin, 
Gansu, China) and Lukuchun (Shanshan, Xinjiang, China) 
with a population from Dzungarian Basin. He found that 
the samples from Dzungarian Basin stood an intermediate 
position between the subspecies E. v. roborowskii 
and nominative form. On the one hand, judging from 
pholidosis, the population from Dzungarian Basin is 
closer to E. v. roborowskii than to the populations from 
Balkhash in Kazakhstan. On the other hand, the dorsal 
and lateral pattern of the population from Dzungarian 
Basin is closer to that of the latter, i.e. nominate form. 
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Interestingly, Szczerbak (1975) referred to the population 
from Dzungarian Basin as nominative form. Our results 
statistically support that the rapid racerunner populations 
in Dzungarian Basin cluster with those from Kazakhstan. 
In fact, biogeographic connections to the Dzungarian 
Basin are to the northwest, with a number of passes 
leading to Kazakhstan. The Dzungarian Gate is the main 
pass through the western ranges, leading to Lake Alaköl 
and Lake Balkhash in Kazakhstan. There are a number 
of taxa that have similar distribution pattern, including 
Phrynocephalus melanurus (Melville et al., 2009), P. 
heliocopus (Solovyeva et al., 2011), and Eremias arguta 
(Szczerbak, 2003). 

Traditionally, the subspecies category is used to 
diagnose geographically distinct populations that were 
thought to be in the early stages of speciation (Mayr, 
1963; Mayr and Ashlock, 1991). Despite some criticisms 
surrounding the concept of subspecies (e.g., Burbrink 
et al., 2000; Patten and Unitt, 2002; Zink, 2004), recent 
studies in which researchers used multiple criteria (e.g., 
morphological, behavioral, and genetic characters) 
have confirmed that many subspecies are evolutionarily 
definable entities (e.g., Gavin et al., 1999; Phillimore 
and Owens, 2006; Guo et al., 2012). In Central Asia, the 
clinal and discrete variation among E. velox populations 
from different regions have been documented (Szczerbak, 
1975; Chirikova, 2004). As noted by Szczerbak (1975), 
the most deviating populations of E. velox were found 
in Fergana Valley and neighborhood of Badkhyz desert, 
reaching subspecies level. Chirikova (2004) studied 
geographic variations of E. velox populations from 
Kazakhstan and Uzbekistan. She found some variations 
in phenotypical traits between populations from the 
two regions. In contrast with the nominate form in 
Kazakhstan, populations from the southern Aral Sea 
region in Uzbekistan have very light elongated spots on 
the middle of the back in adults. Even some individuals 
from the southern Aral Sea region lack this kind of dorsal 
pattern completely. Our results demonstrate that the 
populations in southern Aral Sea region of Uzbekistan 
correspond with distinct, well-supported, and moderately 
divergent mtDNA haplotype clade (clade VI) . The mean 
uncorrected level of sequence divergence between clades 
VI and VII falls within 3.3%. In support of Chirikova 
(2004), we recommend that the populations of rapid 
racerunner from southern Aral Sea region of Uzbekistan 
be assigned to a novel subspecies, which will facilitate an 
effective short-cut for estimating patterns of intraspecific 
genetic diversity of E. velox. 


5. Conclusions 


The mtDNA data confirms the subspecific status for 
three of the four morphology-defined subspecies: 
E. v. roborowskii, E. v. caucasia, E. v. velox (sensu 
stricto). Specifically, we corroborate that there are two 
subspecies occurring in China, i.e., E. v. velox and E. v. 
roborowskii. We recommend a novel subspecies for the 
phenotypically and genetically distinct populations in 
Uzbekistan previously referred to as E. v. velox. We also 
support species-level recognition for each of the three 
independently evolving lineages from Iranian Plateau 
traditionally assigned to Æ. v. velox. We do so recognizing 
that more comprehensive geographic sampling and 
investigation of additional molecular markers are required 
to clarify hypothesized species/subspecies boundaries, 
taxonomic status and evolutionary history of the E. velox 
complex. 
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Appendix A List of analyzed specimens, their geographical origin and the GenBank accession number. 


Accession number 


Taxon Bee Geographical origin Latitude, longitude Sample No. a Seb 12S tRNA Reference 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo475 H1 KF999318 KF999405 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo476 Z1 KF999319 KF999406 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo478 Zi KF999320 KF999407 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo479 Z1 KF999321 KF999408 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo481 Zi KF999321 KF999409 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo483 H2 KF999323 KF999410 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo484 H3 KF999324 KF999411 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E Guo485 Z1 KF999325 KF999412 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E QY1015 Zi KF999396 KF999483 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E QY1016 TUL KF999397 KF999484 This study 
E. v. velox 14 Huocheng, Xinjiang, China 44.02°N, 80.81°E QY1017 Z1 KF999398 KF999485 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo520 Z1 KF999326 KF999413 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo521 Z1 KF999327 KF999414 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo522 Ji KF999328 KF999415 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo523 Z1 KF999329 KF999416 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo524 Zi KF999360 KF999447 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo525 J2 KF999361 KF999448 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo526 J3 KF999362 KF999449 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo527 J4 KF999363 KF999450 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo528 J5 KF999364 KF999451 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo529 Z1 KF999365 KF999452 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo530 TAI KF999366 KF999453 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo531 J6 KF999367 KF999454 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo532 Z1 KF999368 KF999455 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo533 Z1 KF999369 KF999456 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo534 J7 KF999370 KF999457 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo535 Z1 KF999371 KF999458 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo536 J8 KF999388 KF999475 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo537 J9 KF999389 KF999476 This study 
E. v. velox 15 Jinghe, Xinjiang, China 44.55°N, 82.59°E Guo538 Z1 KF999390 KF999477 This study 
E. v. velox 16 Fukang, Xinjiang, China 44.41°N, 87.86°E WGXG08103 Z1 KF999330 KF999417 This study 
E. v. velox 17 Dabancheng, Xinjiang, China 43.34°N, 88.22°E WGXG08525 Zi KF999384 KF999471 This study 
E. v. velox 17 Dabancheng, Xinjiang, China 43.34*N, 88.22°E WGXG08526 wil KF999385 KF999472 This study 
E. v. velox 17 Dabancheng, Xinjiang, China 43.37°N, 88.14°E WGXG08536 Zi KF999386 KF999473 This study 
E. v. velox 17 Dabancheng, Xinjiang, China 43.37*N, 88.14°E WGXG08537 Z1 KF999387 KF999474 This study 
E. v. velox 17 Dabancheng, Xinjiang, China 43.37°N, 88.14°E WGXG08538 W2 KF999359 KF999446 This study 
E. v. velox 17 Dabancheng, Xinjiang, China 43.37°N, 88.14°E WGXG08539 W2 KF999358 KF999445 This study 
E. v. velox 18 Karamay, Xinjiang, China 45.27°N, 85.03°E GP07080433 KL1 KF999404 KF999491 This study 
E. v. velox 19 Qitai , Xinjiang, China 44.29°N, 90.11°E WGXG08116 Ql KF999331 KF999418 This study 
E. v. velox 19 Qitai , Xinjiang, China 44.29°N, 90.11°E WGXG08118 Q2 KF999399 KF999486 This study 
E. v. velox 19 Qitai , Xinjiang, China 44.29°N, 90.11°E WGXG08119 Zl KF999400 KF999487 This study 
E. v. velox 19 Qitai , Xinjiang, China 44.29°N, 90.11°E WGXG08120 Z1 KF999401 KF999488 This study 
E. v. velox 19 Qitai , Xinjiang, China 44.29°N, 90.11°E WGXG08121 Q3 KF999402 KF999489 This study 
E. v. velox 19 Qitai , Xinjiang, China 44.29°N, 90.11°E WGXG08122 Q4 KF999403 KF999490 This study 
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: Accession number 
Taxon Locality Geographical origin Latitude, longitude Sample No. Haplotype Reference 
No. type cyt b 12S rRNA 
E. v. velox 20 Tacheng, Xinjiang, China 46.64°N, 82.95°E JF1293 TAL  KF999395 KF999482 This study 
E. v. velox 21 Fuyun, Xinjiang, China 44.42°N, 88.81°E JF1383 Z1 KF999393 KF999480 This study 
E. v. velox 22 Gaoquan, Xinjiang, China 44.48°N, 84.17°E JF1230 Z1 KF999394 KF999481 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E WGXG08202 Tl KF999332 KF999419 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E WGXG08203 T2 KF999333 KF999420 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E WGXG08204 Tl KF999334 KF999421 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E WGXG08205 T3 KF999335 KF999422 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E WGXG08206 Tl KF999336 KF999423 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E WGXG08207 T4 KF999337 KF999424 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E LJE0707006 T5 KF999372 KF999459 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E LJE0707005 T6 KF999373 KF999460 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E LJE0707004 T7 KF999374 KF999461 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E LJE0707003 T8 KF999375 KF999462 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E LJE0707002 T9 KF999376 KF999463 This study 
E. v. roborowskii 23 Turpan, Xinjiang, China 42.86°N, 89.19°E LJE0707001 T10 KF999377 KF999464 This study 
E. v. roborowskii 24 Jiaohe Ruins, Xinjiang, China 42.95°N, 89.07°E WGXG08208 T11 KF999338 KF999425 This study 
E. v. roborowskii 25 Toksun, Xinjiang, China 42.99°N, 88.75°E WGXG08258 T12 KF999339 KF999426 This study 
E. v. roborowskii 25 Toksun, Xinjiang, China 42.99°N, 88.75°E WGXG08259 T12 KF999379 KF999465 This study 
E. v. roborowskii 25 Toksun, Xinjiang, China 42.99°N, 88.75°E WGXG08260 T12 KF999379 KF999466 This study 
E. v. roborowskii 25 Toksun, Xinjiang, China 42.99°N, 88.75°E WGXG08261 T12 KF999380 KF999467 This study 
E. v. roborowskii 25 Toksun, Xinjiang, China 42.99°N, 88.75°E WGXG08262 T13 KF999381 KF999468 This study 
E. v. roborowskii 25 Toksun, Xinjiang, China 42.99°N, 88.75°E WGXG08263 T12 KF999382 KF999469 This study 
E. v. roborowskii 25 Toksun, Xinjiang, China 42.99°N, 88.75°E WGXG08264 T12 KF999383 KF999470 This study 
E. v. caucasia 26 Baskunchak, Astrakhan, Russia 46.98°N, 48.21°E Guol266 R1 KF999340 KF999427 This study 
E. v. caucasia 27 D River, Astrakhan, 47 gion 46.57*B Gu0l268 RI KF999341 KF999428 This study 
E. v. caucasia 27 Aen River, Astrakhan, 47 gion, 46.57°E Guo1269 RI KF999342 KF999429 This study 
E. v. caucasia 27 ke a River, Astrakhan, 47 gion 46.57°R Guol270 R2  KF999343 KF999430 This study 
E. v. caucasia 27 a a River, Astrakhan, 47 gion 46,57°E Guol271 RI KĶKF999344 KF999431 This study 
E. v. caucasia 27 A a River, Astrakhan, 47 gion. 46,57°E Guol272 RI KF999345  KF999432 This study 
E. v. caucasia 27 = a River, Astrakhan, 47 910N,46.57°E Guol273 RI KF999346 KF999433 This study 
E. v. caucasia 27 eMe a River, Astrakhan, 47 91]oN,46.57°E Guo1274 R3 KF999347 KF999434 This study 
E. v. caucasia 28 be River, Astrakhan, 46 soon, 47.88°E Guol285 R3 KF999350 KF999437 This study 
E. v. caucasia 28 k River, Astrakhan, 46 soon, 47.88°E Guo1286 R6 KF999349 KF999436 This study 
E. v. caucasia 28 bai River, Astrakhan, 46 soon, 47.88°E Guol287 R7 KF999348 KF999435 This study 
E. v. caucasia 28 ka River, Astrakhan, 46 soon, 47.88°E Guol289 RI KF999357 KF999444 This study 
E. v. caucasia 28 Ley River, Astrakhan, 4¢ soon, 47.88°E Guol290 RI KF999356 KF999443 This study 
E. v. caucasia 29 Kalmykia, Russia 46.67°N, 46.29°E Guol276 R4 KF999355 KF999442 This study 
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Taxon Locality Geographical origin Latitude, longitude Sample No. Haplotype eee Reference 

No. type cyt b 12S rRNA 

E. v. caucasia 29 Kalmykia, Russia 46.67°N, 46.29°E Guol277 R1 KF999354 KF999441 This study 

E. v. caucasia 29 Kalmykia, Russia 46.67°N, 46.29°E Guol278 R5 KF999353 KF999440 This study 

E. v. caucasia 29 Kalmykia, Russia 46.67°N, 46.29°E Guol279 Rl KF999352 KF999439 This study 

E. v. caucasia 29 Kalmykia, Russia 46.67°N, 46.29°E Guol284 R1 KF999351 KF999438 This study 

E. v. caucasia 30 Dagestan, Russia 43.61°N, 46.34°E Guol711 R8 KF999391 KF999478 This study 

E. v. caucasia 30 Dagestan, Russia 43.61°N, 46.34°E Guol712 R1 KF999392 KF999479 This study 

Eremias sp. 1 1 North Tehran, Iran — 02-74 El JQ690169 JQ690100 ý 

Eremias sp. 2 3 SE Golestan National park, Iran - Smp221 E3 JQ690175 JQ690106 

Eremias sp. 2 6 S Naishaboor, N Khorasan, Iran — EMP188 E6 JQ690187 JQ690119 

Eremias sp. 3 2 SE Tehran, Iran — ERP266 E2 JQ690172 JQ690131 

Eremias sp. 3 4 10 km S Jajarm town, Iran - Smp268 E4 JQ690171 JQ690102 

Eremias sp. 3 5 15 km SW Sabzevar, Iran - Smp15 E5 JQ690183 JQ690115 

E. velox 7 S Aral Sea, Uzbekistan - Ev2 UZ1 JQ690216 JQ690149 

E. velox 7 S Aral Sea, Uzbekistan - Ev3 UZ2 JQ690217 += JQ690150 

E. velox 7 S Aral Sea, Uzbekistan - Evl UZ3 JQ690218 JQ690151 

E. velox 7 S Aral Sea, Uzbekistan - Ev4 UZ4 JQ690219 JQ690152 . 

E. velox 7 S Aral Sea, Uzbekistan - Ev6 UZ4 = JQ690221 JQ690154 ý 

E. velox 7 S Aral Sea, Uzbekistan - Ev7 UZS JQ690222 JQ690155 j 

E. velox 7 S Aral Sea, Uzbekistan — Ev8 UZ6  JQ690223 JQ690156 * 

E. velox 7 S Aral Sea, Uzbekistan - Ev5 UZ6  JQ690220 JQ690153 

E. velox 7 S Aral Sea, Uzbekistan - Ev9 UZ7 JQ690224 JQ690157 

E. velox 7 S Aral Sea, Uzbekistan - Ev11 UZ8 JQ690225 JQ690158 i 

E. velox 7 S Aral Sea, Uzbekistan - Ev12 UZ9 JQ690226 JQ690159 * 

E. velox 7 S Aral Sea, Uzbekistan - Ev13 UZ10 JQ690227 JQ690160 

E. v. caucasia 9 Extreme W Kazakhstan — D4 KZW1 JQ690233 JQ690166 

E. v. velox 8 W Aral Sea, Kazakhstan - M90 KZ1 JQ690212 + JQ690145 

E. v. velox 10 SW Balkhash Lake, Kazakhstan - K04-6 KZ2 =JQ690215 JQ690148 

E. v. velox 10 SW Balkhash Lake, Kazakhstan - D10 KZ3 JQ690228 JQ690161 

E. v. velox 10 SW Balkhash Lake, Kazakhstan - D9 KZ4  JQ690229 JQ690162 

E. v. velox 10 SW Balkhash Lake, Kazakhstan — RL26 KZ5  JQ690230 JQ690163 

E. v. velox 10 SW Balkhash Lake, Kazakhstan - RL28 KZ6 JQ690231 JQ690164 

E. v. velox 10 SW Balkhash Lake, Kazakhstan - RL65 KZ7 JQ690232 JQ690165 i 

E. v. velox 10 SW Balkhash Lake, Kazakhstan - RL112 KZ8 JQ690234 + JQ690167 . 

E. v. velox 11 NE Aral Sea, Kazakhstan — 04-8 KZ9 JQ690213 JQ690146 * 

E. v. velox 11 NE Aral Sea, Kazakhstan - K04-12 KZ10 JQ690214 JQ690147 . 

E. v. velox 12 Turkistan, S Kazakhstan - T105 KZ11 JQ690206 JQ690139 

E. v. velox 13 E Kazakhstan - M6 KZ12 JQ690204 JQ690137 

E. v. velox 13 E Kazakhstan - U16 KZ13 JQ690205 JQ690138 . 

E. v. velox 13 E Kazakhstan - Pp KZ14  JQ690207 JQ690140 . 

E. v. velox 13 E Kazakhstan - All KZ15  JQ690208 JQ690141 . 

E. v. velox 13 E Kazakhstan - Al KZ16 JQ690209 JQ690142 

E. v. velox 13 E Kazakhstan - U12 KZ17 JQ690210 JQ690143 

E. v. velox 13 E Kazakhstan - M14 KZ18 JQ690211 JQ690144 

E. pesica - Tran - ERP164 - FJ416285 FJ445362 

E. argus - South Korea - LEGO-F415 - JQ086345 JQ086345 i 
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